//
// Pythia6.cxx
//
// Created by TB on 12/6/11.
// Copyright 2011 BNL. All rights reserved.
//

#include <cmath>
#include <cstdlib>
#include <ctime>
#include <iostream>
#include <set>

#include <boost/foreach.hpp>

#include <TPythia6.h>

#include "tppmc/Configuration.h"
#include "tppmc/Enumerations.h"
#include "tppmc/Pythia6.h"

namespace tppmc {
   
   Pythia6::Pythia6() {
   }

   Pythia6::Pythia6(const Configuration& config) {
      Initialise(config);
   }

   Pythia6::~Pythia6() {
   }

   void Pythia6::Initialise(const Configuration& config) {
      TPythia6& EvtGen = *(TPythia6::Instance());

      if(config.PythiaTune() > 0) {
         EvtGen.SetMSTP(5, config.PythiaTune());
         // Do process selection after setting tune parameters
         EvtGen.SetMSEL(0);
         BOOST_FOREACH(int process, config.mProcesses) {
            EvtGen.SetMSUB(process, true);
         } // BOOST_FOREACH
      } // if
      else if(config.pyth_evt_type == QCD2x2){
         // QCD2x2 physics processes
         EvtGen.SetMSEL(0);
         BOOST_FOREACH(int process, config.mProcesses) {
            EvtGen.SetMSUB(process, true);
         } // BOOST_FOREACH
         EvtGen.SetMSTP(81,0);
      } // else if
      else if(config.pyth_evt_type == GENJPSI){
         // Tune from Tevatron (color octet)
         // J/Psi Production
         EvtGen.SetMSEL(61);
         EvtGen.SetMSTJ(11, 3);
         EvtGen.SetMSTJ(22, 2);
         // \todo MSTJ(1) = 0 turns off fragmentation.
         // Make this an options via Configuration.
         EvtGen.SetMSTJ(1, 0);
         EvtGen.SetPARJ(71, 10);
         EvtGen.SetMSTP(2, 2);
         EvtGen.SetMSTP(33, 3);
         EvtGen.SetMSTP(51, 7);
         EvtGen.SetMSTP(81, 1);
         EvtGen.SetMSTP(82, 4);
         EvtGen.SetMSTU(21, 1);
         EvtGen.SetPARP(82, 1.9409);
         EvtGen.SetPARP(89, 200.0);
         EvtGen.SetPARP(83, 0.5);
         EvtGen.SetPARP(84, 0.4);
         EvtGen.SetPARP(90, 0.16);
         EvtGen.SetPARP(67, 2.5);
         EvtGen.SetPARP(62, 1.25);
         EvtGen.SetPARP(64, 0.2);
         EvtGen.SetPARP(85, 1.00);
         EvtGen.SetPARP(86, 1.00);
         EvtGen.SetPARP(88, 0.5);
         EvtGen.SetPARJ(13,0.750);
         EvtGen.SetPARJ(14,0.162);
         EvtGen.SetPARJ(15,0.018);
         EvtGen.SetPARJ(16,0.054);
         EvtGen.SetMSTP(145,0);
         EvtGen.SetMSTP(146,0);
         EvtGen.SetMSTP(147,0);
         EvtGen.SetMSTP(148,1);
         EvtGen.SetPARP(141,1.16);
         EvtGen.SetPARP(142,0.0119);
         EvtGen.SetPARP(143,0.01);
         EvtGen.SetPARP(144,0.01);
         EvtGen.SetPARP(145,0.05);
         EvtGen.SetPARP(146,9.28);
         EvtGen.SetPARP(147,0.15);
         EvtGen.SetPARP(148,0.02);
         EvtGen.SetPARP(149,0.02);
         EvtGen.SetPARP(150,0.085);
         // Decays turned off below
         // This allows ChiC -> J/Psi+gamma but J/Psi is stable  
         EvtGen.SetMDME(858,1,0);
         EvtGen.SetMDME(859,1,0);
         EvtGen.SetMDME(860,1,0);
         EvtGen.SetMDME(1619,1,1);
         EvtGen.SetMDME(1620,1,0);
         EvtGen.SetMDME(1673,1,1);
         EvtGen.SetMDME(1674,1,0);
         EvtGen.SetMDME(979,1,1);
         EvtGen.SetMDME(980,1,0);
      } // else if
      else{
         std::cerr << "ERROR! pyth_evt_type unrecognized!" << std::endl; 
         exit (-1); 
      } // else
      
      // Set the partonic kT, which is treated as a Gaussia distribution.
      // From the PYTHIA 6 manual:
      //   MSTP(91) : (D = 1) primordial kT distribution in hadron.
      //     = 1 : Gaussian, width given in PARP(91), upper cut-off in PARP(93).
      EvtGen.SetMSTP(91,1);
      EvtGen.SetPARP(91, config.primordial_kt);  
      EvtGen.SetPARP(93, 3. * config.primordial_kt); // Change to x5 cutoff?
      
      // Set minimum parton pT
      EvtGen.SetCKIN(3, config.pyth_parton_pt);
      std::cout << " Pythia min PT in hard scattering set to "
      << TPythia6::Instance()->GetCKIN(3) << " GeV" << std::endl;
      
      // Master switch for final state radiation (0=OFF, 1=ON)
      EvtGen.SetMSTP(71, config.order);
      // Master switch for initial state radiation (0=OFF, 1=QCD, 2=QCD+QED)
      if(config.order==0) {
         EvtGen.SetMSTP(61, config.order);
      } // if
      else {
         EvtGen.SetMSTP(61,2);
      } // else

      if(config.pyth_evt_type not_eq GENJPSI) {
         // Turn off hadronization
         EvtGen.SetMSTP(111,0);
         // No decays
         EvtGen.SetMSTP(41,0);
      } // if
      
      // Set random seed
      EvtGen.SetMRPY(1, static_cast<int>(time(NULL))); // random seed number
      
      // Set up the generator based on the machine type
      switch(config.machine_type){
         case(RHIC_200):
            EvtGen.Initialize("CMS","p","p",200.);
            std::cout << " Initialized pythia for p+p at center-of-mass energy = 200GeV " << std::endl; 
            break;
         case(RHIC_500):
            EvtGen.Initialize("CMS","p","p",500.);
            std::cout << " Initialized pythia for p+p at center-of-mass energy = 500GeV " << std::endl; 
            break;
         case(FNAL_E704_PROTON):
            EvtGen.Initialize("CMS","p","p", std::sqrt(2. * 0.938 * 200.));
            std::cout << " Initialized pythia for p+p at center-of-mass energy = " << std::sqrt(2. * 0.938 * 200.)
            << " GeV" << std::endl; 
            break;
         case(FNAL_E704_PBAR):
            EvtGen.Initialize("CMS","pbar","p", std::sqrt(2. * 0.938 * 200.));
            std::cout << " Initialized pythia for pbar+p at center-of-mass energy = " << std::sqrt(2. * 0.938 * 200.)
            << " GeV" << std::endl; 
            break;
         default: 
            std::cerr << "ERROR - Undefined machine_type!" << std::endl; 
            exit(-1); 
      } // switch
   }
} // namespace tppmc
